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Abstract 

In this work we argue that, in the interpretation of the energy dependence of the 
depth of the shower maximum and of the muon content in high energy cosmic 
ray showers (i?>10^^ eV), other variables besides the composition may play an 
important role, in particular those characterising the first (high energy) hadronic 
collisions. The role of the inelasticity, of the nature of the leading particle, and 
of the particle multiplicity are discussed. A consistent interpretation of existing 
data within a string percolation model implemented in a hybrid, one dimensional 
simulation method is given. 
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1 Introduction 



The composition of cosmic rays at high energies is a matter of controversy. 
Direct measurements are possible only up to E ^ 10^^ eV [1]. Above these 
energies, attempts to infer the effective mass number A of the primary par- 
ticle are based on the measured shower variables and are rather indirect [2]. 
Furthermore, data are relatively scarce at such high energies. 

In fluorescence experiments, the longitudinal shower profile is measured and 
the atmospheric depth of the shower maximum Xmax is usually interpreted as 
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an indication of the nature of the primary particle [3]: heavy nuclei initiated 
showers are expected to develop earlier, due to the larger interaction cross- 
section. 

In ground array experiments, the muon component is widely studied as a pos- 
sible handle on composition [3]. The muon content, or more specifically the 
muon density at a given distance from the shower core, p^, is expected to 
increase by roughly 50% from proton to Fe initiated showers, for the same pri- 
mary energy [4] . Incidentally, the muon component of the shower is a powerful 
tool for the validation of hadronic interaction models at such high energies. 

In this work we argue that, in the interpretation of the energy dependence 
of Xmax and p^, other variables besides A may play an important role, in 
particular those characterising the first (high energy) hadronic collisions. The 
role of the inelasticity K, one minus the fraction of momentum carried by 
the fastest (leading) particle, of the nature of the leading particle, and of the 
charged particle multiplicity are discussed. A consistent interpretation of data, 
within a string percolation model implemented in a hybrid one dimensional 
simulation method, is given. 

This paper is organised as follows. In section 2 the experimental status is 
briefly reviewed. In section 3, the possible role of the main variables char- 
acterising the first hadronic coUisions is anahsed and the models and tools 
relevant for the discussion are presented. In section 4, results on the energy 
dependence of X^ax and are derived in the light of string percolation. In 
section 5 some conclusions are drawn. 



2 Experimental status 

Experimental data on the energy dependence of the depth of the shower maxi- 
mum Xmax and of the muon density 600 m away from the shower core p^(600) 
above 10^'' eV, and their comparison with simulations, are summarised in this 
section. 

Fly's Eye/HiRes results in this energy range show an increase in the X^ax vs. 
E slope, which has been interpreted as a change in composition, going from 
more Fe-like to more proton-like showers [5]. It has been argued [6,7] that other 
models predicting a development of the shower deeper in the atmosphere (an 
increase in X^ax) could also explain this effect. 

Concerning the muon component of the shower, experimental data are rela- 
tively scarce at high energies. Extensive air showers are largely electromagnetic 
dominated, and the separate detection of muons in ground array experiments 
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requires the installation of specific shielded detector units. This has been done 
in KASCADE [8] in the energy range from 10^^ eV to almost 10^'' cV and in 
AGASA above 10^^ eV [9]. Hybrid data from HiRes-MIA has also been pre- 
sented [10]. 

At AGASA, the lateral distribution function of muons above 0.5 GeV was mea- 
sured and combined with the Akeno 1 Km^ array data (threshold 1 GeV) [9]. 
The evolution of Pn{600) with the particle density 5*0(600) (essentially propor- 
tional to the primary energy) was studied and compared with CORSIKA [11] 
simulations using different hadronic models in [4]. It was observed that the 
slope of p/i(600) vs. E in data is fiatter than in simulation with any hadronic 
model and primary composition. Taking cx E^, data gives f3 = 0.84 ± 
0.02 [4,9], while from simulation (3 0.9 {f3 = 0.92 (0.89) for protons with 



QGSJET (SIBYLL), while for Fe /3 = 0.88 (0.87) [4]). It has been argued [12] 



that these results can be interpreted as a change in composition from heavy 
(at around lO^''''^ cV) to light (at around 10"^^ eV), agreeing with the Fly's Eye 
indication. It has been pointed out in [4] that this interpretation should be 
studied in a wider energy range, as it seems to lead to a composition heavier 
than Fe at lower energies. The HiRes-Mia experiment consisted of the pro- 
totype high resolution Fly's Eye detector and the muon array MIA. Results 
are presented in [10]. The muon data, yielding /3 = 0.73, seem to indicate a 
composition heavier than Fe at 10^^ eV. 



3 Construction of a simple model 

3.1 The composition interpretation 

In order to make the argument simple, we shall use the original Heitler idea [13]: 
the location of the shower maximum is, on the average, related to log E, X^ax ~ 
log E, and the number of charged pions or muons is proportional to some power 
of E, A^TT+TT- ^ ^ E'^, < /3<1. In the case of a nucleus with A nucleons 
colliding in a hadronic collision we shall write. 



where Xi is the average depth of the first collision and Xq is the radiation 
length, and 



X, 



max — 



X, + Xolog{E/A), 



(1) 



N, ~ A{E/A)^, 



(2) 
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or 



logTV^ ~ (1 -/3)logA + /31ogE 



(3) 



We further have, for the elongation rate, 



dXr, 



dlogE 



1 - 



dlog^ 
dlogE 



(4) 



and for the log dependence on E, 
rflogiV,^^^_^^rflogA^^_ 



diog£; 



dlogE 



(5) 



As experimentally [14,5,10] , above 10^^ eV dXmax/ d log E is larger and d log N^/ 
dlogE is slightly smaller, in comparison with lower energies, the conclusion 
is: 



dlogA 
dlogE 



<0, 



(6) 



i.e., for £'>10^^ eV the average mass number A should, in this interpretation, 
decrease with energy. 



3.2 The first hadronic collisions interpretation 



In this interpretation a key role is given to the main variables characterising 
the first hadronic collisions, which are schematically represented in Fig. 1. 
The inelasticity K is the fraction of energy distributed among the produced 
particles (mainly pions) . A fraction of this energy is expected to go into neutral 
pions, which promptly decay to photons. < n > is the average (non-leading) 
multiplicity at the collision energy. The quantity (l—K) represents the fraction 
of energy concentrated in the leading particle, in general assumed to be a 
proton. Pq is the probability of having a leading 7r°. 

As pointed out in [6,15], the effect of changing A is equivalent to the effect of 
changing the average inelasticity K. In this spirit, the fastest particle, carrying 
an energy {1 — K)E, will originate the shower branches that go deeper in the 
atmosphere, and we write 

Xmax ^X,+Xo l0g[(l -K)^], (7) 
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instead of (1), where Xi is the average depth of the first colhsion, Xq is the 
radiation length and Eq a low energy threshold. 




<n> 



Fig. 1. Schematic representation of the first hadronic collisions. A large fraction 
of the energy (1 — K) is taken by the leading particle. Most particles (essentially 
pions, < n > on average) come from the remaining energy, K. The probability Pq 
characterises the type of leading particle. 

Regarding the muon content of the shower, a possible assumption is to say 
that while the energy flows in the (1 — K) direction (see Fig. 1) the number 
of particles (of muons) flows in the K direction. We thus have 

^ oc KE . (8) 



instead of (2). The validity of this assumption will be discussed later. We 
further have. 



dX„ 



dlogE 



d\og{l-K) 
dlogE 



(9) 



instead of (4), and 

dlogN^ _ dlogK 
dlogE d log E 



+ 1 



(10) 



instead of (5). The condition (6) becomes now 
dlogK 



dlogE 



<0 , 



i.e., the inelasticity K has to decrease with the energy. This scenario has 
been explored in [7], where a string percolation model predicting the required 
behaviour of the inelasticity at high energies was discussed. As it will be seen 
below, string percolation affects the energy dependence of K, Pq and < n >. 
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3.3 Matthews-Heitler toy model 



In [16] a very interesting semi-empirical model for the development of the 
hadronic component of air showers in analogy with Heitler's splitting approx- 
imation of electromagnetic cascades [13] was presented. Not considering the 
elasticity, the model gives constant and reasonable values for dX^ax/dlogE 
and d log iV^/ d log E. Taking into account that in hadronic interactions a sig- 
nificant fraction of the total energy may be carried away by a single "leading" 
particle (i.e., 1 — K > 0), the model predicts for a slope decreasing with 
K - exactly the opposite of (8). It is interesting to discuss this discrepancy. 

In [16] the essential point is that the production of 7r°'s, which decay into 
photons before participating in the cascade, strongly decreases the muon pro- 
duction - about 1/3 of the energy KE is lost for muon production, as it goes 
into the electromagnetic component of the shower. On the other hand, the 
leading particle, assumed in general to be a proton, goes on producing par- 
ticles, thus increasing the number of muons. However, if tt^'s are themselves 
leading particles - Pq > 0, as it happens in percolation models - the argu- 
ment is no longer valid. In fact, when estimating muon production both the 
inelasticity K and the probability Pq of having a 7r° as a leading particle are 
relevant. This will be seen next. 



3.4 Generalised Matthews-Heitler model 

The model of [16] was implemented and generalised. This is not a substitute 
for a detailed simulation, but useful to illustrate the physics involved. Follow- 
ing [16], hadrons interact after traversing one layer of fixed thickness (related 
to the interaction length, Xj ~ 120 g cm^ [16,17]) producing Nch charged 
pions and Nch/'^ neutral pions. A 7r° immediately decays into photons, mean- 
ing that 1/3 of the energy is lost to the electromagnetic shower component. 
Charged pions continue through another layer and interact, until a critical 
energy is reached. They are then assumed to decay, yielding muons, and thus 
— N^± . The version of the model which takes into account the inelasticity 
was implemented. In each branching, the energy was equally shared amongst 
all (non-leading) particles. Nch = 10 was used. 

We modified this model in two ways. Firstly, a decay probability for charged pi- 
ons in each interaction step was explicitly included, replacing the large value of 

the critical energy (20 GeV) used in the original version. This gives a smoother 
and more realistic transition from charged pions to muons. A critical energy of 
1 GeV was then used. Further, the model was generalised to include the prob- 
abihty Pq of having a leading 7r° in a given collision, with an energy (l-K)E 
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lost into the electromagnetic branch in that collision. The obtained number of 
muons as a fimction of the inelasticity K (for E — 10^^ eV) is shown in Fig. 2, 
for different values of Pq. 

As observed, in this toy model the slope depends critically on the nature 
of the leading particle: changing the probability Pq of a leading neutral pion 
from to 1/3 inverts this slope. If the inelasticity is relatively low, even a 
moderate fraction of leading neutral pions may cause a decrease of the number 
of produced muons. Thus, one easily moves from the behaviour predicted in 
the Matthews-Heitler model (section 3.3) to the first hadronic collisions model 
(section 3.2, equation (8)). It should however be noted that in this simplified 
approach we do not take into account the fact that the energy spectrum of the 
leading particle will depend on its nature: leading neutral pions may be softer 
than leading protons. More detailed simulations will be presented below. 
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Fig. 2. The toy model prediction for the the number of muons as a function of the 

inelasticity is shown, for different values of the probability Pq of having a leading 
7r° and E = 10^^ eV. The lines are only to guide the eye. 



3. 5 Hybrid shower simulations 



In order to validate the argument and derive predictions for the behaviour 

of Xmax and as a function of the energy, we need a shower simulation 
tool which is detailed enough to produce reliable results, and fast and flexible 
enough to allow the ad — hoc introduction of the predictions of the percolation 
model that will be used. 
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The hybrid, one dimensional simulation method described in [18] was used. It 
is a fast, one dimensional calculation, which provides predictions for the total 
number of charged particles and muons above several energy thresholds along 
the shower axis, as well as for the fluctuations of the electromagnetic and 
hadronic components of the shower. The method is based on precalculated 
showers (sub-threshold particles are treated with a library of profiles based 
on pre-simulated pion-initiatcd showers) and a bootstrap procedure to extend 
the shower library to high energy. The SIBYLL hadron interactions model was 
used. 

In Fig. 3, the number of muons {E > 1 GeV) as a function of K is shown, for 
E = 10^^ eV, and for different values of the fraction of leading vr^'s, Pq- This 
confirms the indication obtained above with the toy model: The slope of 
vs. K depends strongly on Pq. 

3000 



2800 



2600 
^ 2400 
2200 
2000 



' 0.5 0.55 0.6 0.65 0.7 0.75 0.8 

K 

Fig. 3. The hybrid simulation prediction for the the number of muons (E > 1 GeV) 
as a function of the inelasticity is shown, for different values of the probability Pq 
of having a leading 7r° and E = 10^^ eV. The lines are only to guide the eye. 

Studies performed with both QGSJET and SIBYLL show that the prediction 
of CORSIKA for vs. K is essentially flat. In QGSJET, the fraction of 
leading tt^'s in the flrst hadronic collision goes roughly from 10% to 15% 
as the primary energy goes from 10^^ eV to 10^° eV, but their spectrum is 
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relatively soft. The variation of Pq in SIBYLL is very small (of the order of 
1% for the same energy range). 

The combined effect of K and Pq has thus an effect on muon production, and 
this is particularly true if the inelasticity has relatively low values. We now 
need a consistent model with predictions for the values of these variables. 

3. 6 String percolation model 

Essentially, all existing high energy strong interaction models based on QCD, 
and QCD evolution, predict an increase with energy - not a decrease - of the 
inelasticity K [19]. The same is true for the hadronic generators SIBYLL [20] 
and QGSJET [21], used in cosmic ray cascade analysis. This happens because 
evolution in the energy implies transfer of energy from valence partons or 
strings, or from bare Pomeron diagrams, to sea partons or strings, or to multi- 
Pomeron contributions. 

However, in models with percolation of partons or strings, one expects the 
inelasticity K to decrease with energy above the percolation threshold [7]. In 
the framework of the Dual String Model [22] - but we beheve the argument is 
more general - what happens at low energy is the transfer of energy from the 
valence strings to sea strings (and K increases), while at higher energy the 
strings start to overlap and a cumulative effect occurs: the length in rapidity 
of fused (percolated) strings is larger [23]. At some stage, close to percola- 
tion threshold, the percolating strings take over the valence strings, and from 
then on K decreases with the energy. Percolation is, in fact, a mechanism for 
generating fast leading particles. 

In Fig. 4 we show the energy dependence of K in the case of our string percola- 
tion model [7], in comparison with K determined from QCD inspired models, 
without percolation (see for instance [24] ). As discussed in [7], this model 
was tuned taking into account multiplicity data from accelerator energies and 
places the percolation threshold at £" ~ 10^^ eV. In the reference, predictions 
for the inelasticity and multiplicity behaviour were made and their effect on 
Xmax data was discussed. Above the percolation threshold the inverted slope 
of K vs. E is clearly visible. In the relevant energy range, K varies roughly 
from 0.8 to 0.55. Percolation will also affect the particle multiplicity, and we 
follow the treatment in [7], where we see that this effect can be introduced 
through the colour summation reduction factor F{r]), where r] is the transverse 
density. 

Let us now turn to the evolution with energy of the probability Pq of having a 
non-baryonic primary. As energy increases, the number of sea strings increases. 
These are, however, low energy strings, and, without percolation, the initial 
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Fig. 4. Energy dependence of the average inelasticity K. Conventional QCD model 
results are from [24]: dotted curve is SIBYLL simulation; dashed-dotted curve is 
QGSJET simulations; dashed is a model of [24]; full line is our percolation model. 

proton remains dominant. When sea strings percolate, larger "sea clusters" 
are formed, and other particles can be produced. It is a prediction of our 
percolation model that above the percolation threshold, because sea strings 
are of the type quark-antiquark, the probability of having a leading 7r° will 
tend to 1/3. 



4 Results 



We now have a string percolation model which predicts modifications on the 
variables K, Pq and < n > characterising hadronic collisions above a percola- 
tion threshold. Inserting these predictions into the hybrid shower simulation 
tool described in section 3.5, results for the behaviour of N/^ and Xmax as a 
function of the energy can now be derived and compared with the available 
data (see section 2). 

The percolation prediction for the inelasticity K was introduced in the hy- 
brid simulations for hadronic coUisions above 10^^ eV. Above the percolation 
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threshold (E ~ 10^^ eV), the probabihty of having a leading 7r° of Pq was set 
to 1/3 and the percolation multiplicity reduction factor F{r]) was introduced. 

In Fig. 5 our results for A^^ vs. E are shown and compared with Akeno/AGASA 
data [9,25]. A value oi P — 0.83 is obtained from simulation, in very good 

agreement with the measured value of f3 = 0.84 ± 0.02. In this study a pri- 
mary proton was considered. The result is however more general, as it has 
been shown in [4] that, for fixed composition, the value of /3 is essentially the 
same for proton and for Fe. The curve shown in the figure corresponds to 
shifting our proton result, using equation (2), to an intermediate and constant 
composition A ~ 20. The result of f3 = 0.73 quoted by HiRes-Mia can hardly 
be accommodated in this model, as it would imply values of Pq going well 
beyond the expectations. 

It is interesting to note that at very high energies {E>10^^ eV) the 7r°'s starting 
to participate in the shower (i.e. to interact) and hence contribute to the 
production of muons. This causes an increase in the slope. A slight effect is 
already seen in our simulations. Finally, as muon data is often presented as 
muon density p^(600) vs. E, it is worth noting that CORSIKA [11] was used 
to check that p^(600) behaves much like A*"^ as a function of the energy (with 
A^^ slightly steeper, but differences of the order of 0.02). 

In Fig. 6 the results obtained for X^ax vs. E are shown and compared with 
existing data. A composition very similar to the one above was used. We see 
that the obtained curve is reasonably consistent with data and with the result 
presented in [7] using a simple model similar to the one described in section 3.2. 

We thus conclude that a percolation model with predictions for the evolution 
with energy of the inelasticity K, the multiplicity < n > and the fraction Pq 
of leading neutral pions, in high energy hadronic collisions, can contribute to 
explain consistently the behaviour observed in data, both for A^^ vs. E and 

Xmax VS. E. 



5 Summary and conclusions 

In this paper we have shown that, in the interpretation of the energy depen- 
dence of the depth of the shower maximum and of the muon content in high 
energy cosmic ray showers (£'>10^^ eV), other variables besides the compo- 
sition may play an important role, in particular those characterising the first 
(high energy) hadronic collisions. The inelasticity K, the nature of the leading 
particle Pq, and the particle multiplicity < n > were discussed. 

We developed a model which includes a hybrid one-dimensional Monte Carlo 
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Fig. 5. The number of muons as a function of the primary energy. The results of 
our model (full line) are shown and compared with Akeno/AGASA data (symbols). 

simulation [18] and expectations from the percolation scenario [7]. For the 
first time the predictions of a percolation model for the different variables 
were included systematically and used to derive predictions on the behaviour 
of the depth of the shower maximum Xmax and on the number of muons in 
the shower A^^ as a function of the energy. We were able to explain the data 
in a consistent way. The model reasonably describes the trends seen in data. 

It shoud be noted however that the model (in particular the behaviour of the 
number of strings as a function of the energy [7]) is tuned for high energy, the 
region where the percolation threshold arises (£'>10^^ eV). The next step will 
be to create a full Monte Carlo simulation including percolation, adapted to 
the low and high energy regions. This is particularly important for the case 
of N^, due to the influence of Pq- A detailed Monte Carlo calculation is thus 
required for a better understanding of data at all energies. This will allow 
to further test the predictions of percolation, namely the relative decrease of 
multiplicity, the non monotonical behaviour of the inelasticity and the fast 
growth of Pq with energy, around 10^^ eV. 

Finally, and concerning the leading tt^'s, it is important to realise that above 
10^^ eV they become active in the cascade, thus generating more hadrons. The 
(3 parameter is expected to increase, and this is already seen in our simulations. 
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Fig. 6. The depth of the shower maximum as a function of the primary energy. The 
results of our model (full line) are shown and compared with HiRes data (symbols). 

However, a proper treatment of this problem requires, once again, a detailed 
Monte Carlo on percolation. More data at extremely high energies are also 
needed. 
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